import commands
from Parameters import *
from scipy.stats import mstats
import numpy.matlib as np

class R2:
	def __init__(self, parameters):
		self.parameters = parameters
		self.setNull()
	
	def runMSQuick_Null( self, outFile='/tmp/ms.out' ):
		newCommand = '~/bin/msQuick' + self.parameters.getMSString_Null() + ' -seeds 1 1 1 | cut -f 1'
		print '\nRunning\n', newCommand
		msOut = commands.getoutput( newCommand )
		
		f = open( outFile, 'w')
		f.write( msOut )
		
		return msOut

	def runMSQuick_Alternative( self, altDemography, outFile='/tmp/ms.out' ):
		newCommand = '~/bin/msQuick' + self.parameters.getMSString_Alternative( altDemography ) + ' -seeds 1 1 1 | cut -f 1'
		print '\nRunning\n', newCommand
		msOut = commands.getoutput( newCommand )
		
		f = open( outFile, 'w')
		f.write( msOut )
		
		return msOut

	def setNull(self):
		msOut = self.runMSQuick_Null()
		
		self.myValues = []
		firstLine = True
		
		for line in msOut.split('\n'):
			if firstLine:
				firstLine = False
				continue

			v = float(line)
			self.myValues.append(v)

		self.getRejectionValues(0.05)
		
	def getRejectionValues(self, alpha):
		rValues = mstats.mquantiles( self.myValues, [alpha/2.0, 1 - alpha/2.0])
		
		self.myRejectLower = rValues[0]
		self.myRejectUpper = rValues[1]
		return rValues

	def getPower(self, altDemography ):
		msOut = self.runMSQuick_Alternative(altDemography)
		
		altValues = []
		firstLine = True
		
		for line in msOut.split('\n'):
			if firstLine:
				firstLine = False
				continue

			v = float(line)
			altValues.append(v)

		myRejected = self.testAltHypotheses( altValues )
		
		numRejected = 1.0 * np.sum( myRejected )
		power = 1.0 * numRejected / self.parameters.myPopulation.numPowerSimulations
		beta = 1.0 - power
		
		return power
	
	def testAltHypotheses(self, inList):
		rValue = None
		rValue = np.zeros( len(inList), int )
		
		for i, value in enumerate(inList):
			if self.isRejected(value):
				rValue[0,i] = 1
				
		return np.ravel(rValue)
		
	
	def isRejected(self, value):
		if (value < self.myRejectLower) or (value > self.myRejectUpper):
			return True
		return False


#p = Population(NUM_OF_CHROMOSOMES, N, LENGTH, MU, 0.0, S)
#d = Demography(N)
#simParams = SimulationParameters( p, d )
#print simParams
#
#r = R2()
#temp = r.setPower(simParams)
#
#altDemography = Demography(N)
#altDemography.addEpoch(t, NEW_SIZE, 0.0)
#temp = r.getPower(altDemography)
#print temp
